input data



source(file = '../inputProcessing/date_input.r')

rep <- 1e4
rep_sim <- 1e3

for (date_num in 1:length(dates_input)){
  # date_num = 1
  
  date_week_finishing <-  as.Date(dates_input[date_num],format = '%Y-%m-%d')
  output <- list()  
  output2 <- list()
  
  for (Mdata in c('A','G')){
    # Mdata = 'A'
    
    inputs<- readRDS(file=paste0(
      '../../Mobility.2.0.Save/Saved_file/inputProcessing/Rdata/',
      date_week_finishing,'_inputs_from_',Mdata,'.rds'))
    
    for (si in 1:2){
      # si=1
      
      ########################################
      ####inputs
      D <- inputs$D
      M <- inputs$M
      Ot <- inputs$Ot[[si]]$Ot
      W <- inputs$W[[si]]$W
      H <- inputs$H
      SI <- inputs$SI[[si]]$SI
      delta_id <- inputs$delta_id
      #
      country <- names(D)[-1]
      
      mD <- as.matrix(D[,-1])
      N_geo <- ncol(D)-1
      N_days <- nrow(D)
      
      N_week <- nrow(D)/7
      
      # prerp MCMC
      # epi
      R0 <- 5
      # for risk
      b <- 0
      o <- 1
      
      theta0 <- c(rep(R0,N_geo),
                  rep(b,N_geo),
                  rep(o,N_geo))
      
      prior_theta <- matrix(c(rep(c(0,5),N_geo),
                              rep(c(-1e2,1e2),N_geo),
                              rep(c(0,1e5),N_geo)),
                            length(theta0),2, byrow=TRUE)
      
      # parameter names
      f0 <- function(x) paste0('R0_',x)
      f1 <- function(x) paste0('beta_',x)
      f2 <- function(x) paste0('Over_',x)
      n_t<- c(sapply(country,f0), sapply(country,f1), sapply(country,f2))
      
      # sd dev for proposal
      sigma <- rep(1e-1,length(theta0))
      
      
      ########################################
      # useful functions
      # useful functions
      sapply(paste0('Rscript/std/',(list.files('Rscript/std/'))),FUN = source)
      
      ########################################
      #run MCMC
      
      #check
      # res <- MCMC_iter(iter = rep, theta0 = theta0, s = sigma)
      res <- MCMC_full(iter = rep, 
                       theta0 = theta0,
                       s = sigma, 
                       repli_adapt = 10, 
                       within_iter = rep/10)
      
      ########################################
      ## check convergence
      
      print('########################################')
      print(paste0('check convergence',date_week_finishing,' - ',Mdata,' - ',si)) 
      
      
      Acc <- colSums(diff(res$theta)!=0)/(rep-1)
      if(rep >1e2) f <- round(seq(1,rep,length.out = 1e2))
      plot(res$logL[,1],main=paste0('DIC ',res$DIC[1],', P ',res$DIC[2]))
      layout(matrix(1:4,2,2))
      for (i in 1:length(theta0)){
        plot(res$theta[,i],
             main = paste0(n_t[i],' - ',round(Acc[i]*100)))#,       ylim = prior_theta[i,])
      }
      
      # apply(res$theta,2,quantile,c(.5,.025,.975))
      
      
      ncountry <- country
      ncountry[which(ncountry %in% 'United_Kingdom')] <- 'UK'
      ncountry[which(ncountry %in% 'United_States_of_America')] <- 'USA'
      
      R0s <- apply(res$theta[,1:N_geo],2,quantile,c(.5,.025,.975))
      betas <- apply(res$theta[,N_geo+(1:N_geo)],2,quantile,c(.5,.025,.975))
      over <- apply(res$theta[,2*N_geo+(1:N_geo)],2,quantile,c(.5,.025,.975))
      
      ## margin for side 2 is 7 lines in siz
      layout(1)
      par(mar = c(7,4,4,2) + 0.1) ## default is c(5,4,4,2) + 0.1
      
      errbar(1:N_geo,R0s[1,],R0s[2,],R0s[3,],
             xlab = '', ylab = TeX('R_0'), bty = 'n',xaxt = "n",ylim=c(0,10))
      # xlab = '', ylab = 'R',ylim = c(0,3), bty = 'n',xaxt = "n")
      
      lines(c(1,N_geo),rep(1,2), col = 'red')
      
      axis(1, at=1:N_geo, labels=ncountry,las=2)
      
      ####
      errbar(1:N_geo,betas[1,],betas[2,],betas[3,],
             xlab = '', ylab = TeX('$\\beta$'), bty = 'n',xaxt = "n")
      # xlab = '', ylab = 'R',ylim = c(0,3), bty = 'n',xaxt = "n")
      
      # lines(c(1,N_geo),rep(1,2), col = 'red')
      
      axis(1, at=1:N_geo, labels=ncountry,las=2)
      
      
      output[[si]] <- list(resMCMC = res,
                           R0s = R0s, 
                           betas = betas, 
                           over = over
      )
      
      #################################### 
      ## Rt daily and Rt_D2
      ### observed mobility
      
      res_m <- output[[si]]
      rep <- nrow(res_m$resMCMC$theta)
      Rt_daily <- array(NA,dim = c(N_days,N_geo,rep))
      Rt_D2 <- Rt_daily
      m_eff <- Rt_daily
      
      
      f_m_eff <- function(){
        1 + 1/B * log( inputs$H %*% exp( - B * (1-inputs$M) ) )
      }
      
      
      for (j in 1:rep){
        R_daily <- Rt_fun(theta = res_m$resMCMC$theta[j,], x = inputs$M )
        Rt_daily[,,j] <- R_daily
        Rt_D2[,,j] <- inputs$H %*% R_daily
        
        B <- rep(1,nrow(inputs$M)) %*% t(res_m$resMCMC$theta[j,(2-1)*N_geo+ (1:N_geo)])
        m_eff[,,j] <- f_m_eff() # 1 + 1/B * log( inputs$H %*% exp( - B * (1-inputs$M) ) )
      }
      
      ###################################
      ### assumed mobility
      
      
      n_d <- 1e3
      Rt_assumed_mob <- array(NA,dim = c(n_d,N_geo,rep))
      x <-  matrix(seq(0,1,length.out = n_d),n_d,N_geo,byrow = FALSE)
      
      for (j in 1:rep){
        R_daily <- Rt_fun(theta = res_m$resMCMC$theta[j,], x = 1-x )
        Rt_assumed_mob[,,j] <- R_daily
      }
      
      
      
      #####################################
      ## summaries
      
      ### Summary for Rt daily and Rt_D2
      
      
      
      # format outputs
      temp <- D
      temp[,-1] <- NA
      
      results_full_Rt_daily <- list(median_R = temp,
                                    low_R = temp,
                                    up_R = temp,
                                    M = temp)
      
      results_full_Rt_D2 <- list(median_R = temp,
                                 low_R = temp,
                                 up_R = temp,
                                 M_D = temp)
      
      results_meff <- list(median_meff = temp,
                           low_meff = temp,
                           up_meff = temp)
      
      for (i in 1:N_geo){
        
        temp <- apply(Rt_daily[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
        results_full_Rt_daily$median_R[,i+1] <- temp[1,]
        results_full_Rt_daily$low_R[,i+1] <- temp[2,]
        results_full_Rt_daily$up_R[,i+1] <- temp[3,]
        results_full_Rt_daily$M[,i+1] <- M[,i]
        
        temp <- apply(Rt_D2[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
        results_full_Rt_D2$median_R[,i+1] <- temp[1,]
        results_full_Rt_D2$low_R[,i+1] <- temp[2,]
        results_full_Rt_D2$up_R[,i+1] <- temp[3,]
        # results_full_Rt_D2$M_D[,i+1] <- M_D[,i]
        
        temp <- apply(m_eff[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
        results_meff$median_meff[,i+1] <- temp[1,]
        results_meff$low_meff[,i+1] <- temp[2,]
        results_meff$up_meff[,i+1] <- temp[3,]
        
        
      }
      
      
      #############################################
      ### Summary for Rt assumed mobility
      
      
      
      # format outputs
      temp <- data.frame(matrix(NA,nrow = n_d, ncol = 1+N_geo))
      temp[,1] <- x[,1]
      names(temp) <- c('mobility', country)
      
      results_full_Rt_assumed_mob <- list(median_R = temp,
                                          low_R = temp,
                                          up_R = temp)
      
      Thresholds <- output[[1]]$R0s
      f <- c()
      for (i in 1:N_geo){
        
        temp <- apply(Rt_assumed_mob[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
        results_full_Rt_assumed_mob$median_R[,i+1] <- temp[1,]
        results_full_Rt_assumed_mob$low_R[,i+1] <- temp[2,]
        results_full_Rt_assumed_mob$up_R[,i+1] <- temp[3,]
        
        f[1] <- which(results_full_Rt_assumed_mob$median_R[,i+1] < 1)[1]
        f[2] <- which(results_full_Rt_assumed_mob$low_R[,i+1] < 1)[1]
        f[3] <- which(results_full_Rt_assumed_mob$up_R[,i+1] < 1)[1]
        for(j in 1:3){
          if (is.na(f[j])){
            Thresholds[j,i] <- NA
          }else{
            Thresholds[j,i] <-mean(results_full_Rt_assumed_mob$median_R$mobility[f[j]:(f[j]+1)])
          }
        }
      }
     output[[si]]$Thresholds <- Thresholds
      
      #############################################
      # forecast and predictions
      
      ## short-term
      n_pred <- 7
      # choose future mobility
      M_pred <- rbind(M,matrix(M[nrow(M),],n_pred,N_geo,byrow = TRUE))
      # useful functions
      sapply(paste0('Rscript/proj_f/',(list.files('Rscript/proj_f'))),FUN = source)
      res_projection_ST <- mobility_prediction_over2(n_pred = n_pred, rep_sim = rep_sim)
      
      ## summary
      
      temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
                         matrix(NA,n_pred,N_geo))
      names(temp) <- c('dates',country)
      summary_proj <- list(median = temp,
                           low = temp,
                           up = temp)
      for (i in 1:N_geo){
        temp <- apply(res_projection_ST$D_pred[,i,],1,quantile,c(.5,.025,.975))
        summary_proj$median[,i+1] <- temp[1,]
        summary_proj$low[,i+1] <- temp[2,]
        summary_proj$up[,i+1] <- temp[3,]
      }
      
      
      output2[[si]] <- list(results_full_Rt_daily = results_full_Rt_daily,
                            results_full_Rt_D2 = results_full_Rt_D2,
                            results_meff = results_meff,
                            results_full_Rt_assumed_mob = results_full_Rt_assumed_mob,
                            summary_proj_ST = summary_proj)
      
    }
    
    saveRDS(object = output,
            file = paste0('../../Mobility.2.0.Save/Saved_file/FitModels/R0_beta_Over/Rdata/',
                          date_week_finishing,'_output_over_form_',Mdata,'.rds' ))
    saveRDS(object = output2,
            file = paste0('../../Mobility.2.0.Save/Saved_file/FitModels/R0_beta_Over/Rdata/',
                          date_week_finishing,'_output2_over_form_',Mdata,'.rds' ))
  }
}
## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 1"

## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 2"

## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 1"

## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 2"

## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 1"

## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 2"

## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 1"

## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 2"